Predator-prey cycles from resonant amplification of demographic stochasticity 
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In this paper we present the simplest individual level model of predator-prey dynamics and show, 
via direct calculation, that it exhibits cycling behavior. The deterministic analogue of our model, 
recovered when the number of individuals is infinitely large, is the Volterra system (with density- 
dependent prey reproduction) which is well-known to fail to predict cycles. This difference in 
behavior can be traced to a resonant amplification of demographic fluctuations which disappears 
only when the number of individuals is strictly infinite. Our results indicate that additional biological 
mechanisms, such as predator satiation, may not be necessary to explain observed predator-prey 
cycles in real (finite) populations. 
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Predator-prey cycles are one of the most striking phe- 
nomena observed in population biology, and as such, in- 
spire intense discussion among ecologists 0, Q- Cycles 
are also seen in a wide variety of other "host-natural en- 
emy" systems, such as host-pathogen systems - one 
of the most well known examples is measles, epidemics of 
which have been studied for many years In this pa- 
per, we will be concerned with modeling the phenomenon 
of cycles, and will focus on predator-prey systems, for 
concreteness, but our main results will have direct appli- 
cability to other host-natural enemy systems, since they 
can be modeled in a similar way, often using identical 
equations. We also believe that, since the phenomenon 
we describe is quite generic in certain classes of stochastic 
systems, it should be found outside population dynam- 
ics. It seems that the precise mechanism underlying the 
existence of the cycles has not so far been elucidated be- 
cause it involves the analysis of stochastic systems with 
a large, but finite, number of constituents, and also be- 
cause it involves concepts such as resonance, which are 
more familiar to physicists than biologists. 

Among the numerous hypotheses put forward to ex- 
plain cycles, perhaps the simplest is that cycles arise 
directly from predator-prey interactions. Within this 
conceptual framework, theoretical modeling of cycles 
has traditionally been developed using deterministic 
population-level models (PLMs). Such discussions begin 
with Volterra-like equations, which are coupled differen- 
tial equations for the predator and prey densities. Equa- 
tions of this type encapsulate the simplest processes of 
predator and prey mortality, prey reproduction and com- 
petition, and predation. Surprisingly these models do not 
predict stable cycles: additional biological mechanisms, 
such as predator satiation, need to be included within 
the framework of differential equations to give cycles . 
It seems puzzling that cycles, which are so easy to under- 
stand intuitively, can only be described mathematically 
in models which include these more subtle mechanisms. 



In order to probe this issue, we shall take a different ap- 
proach here, and describe the predator-prey system using 
an individual level model (ILM). The individuals, which 
are either predators or prey, are acted upon by simple 
stochastic processes of mortality, reproduction, and pre- 
dation. We are able to derive an exact description of 
this model when the number of individuals is large and 
finite. We find that the predator and prey numbers un- 
dergo large cycles, just as one would expect intuitively. 
The cycles, which arise from a novel resonance effect, dis- 
appear only when the number of individuals is taken to 
be strictly infinite, that is, when the PLM is recovered. 
From a statistical physics viewpoint, we would term the 
PLM a mean field theory of the underlying "microscopic" 
ILM which includes statistical fluctuations. 

Predator-prey cycles observed in nature will have a 
stochastic component — this will affect both their am- 
plitude and phase. Therefore care must be taken in aver- 
aging over replicates. A direct average of the population 
densities from different replicates will result in a constant 
average density since, in the absence of an external "forc- 
ing" , there is a lack of synchrony between the cycles from 
different replicates. This fact is crucial when modeling 
predator-prey cycles. A given PLM is written in terms 
of a population density, which can be thought of as the 
result of an average of the population numbers from a 
large number of ILM replicates. If a given ILM shows 
oscillatory behavior, such cycles will be lost in the mod- 
eling transition to a PLM. Thus, it is necessary to study 
quantities such as the autocorrelation function and power 
spectrum arising from an ensemble of ILMs, in order to 
determine the presence and properties of predator-prey 
cycles. 

The specific ILM we study in this paper is a non-spatial 
stochastic model. At a given time, a realization of the 
ILM consists of n individuals of species A (the preda- 
tors) and m individuals of species B (the prey) . Since we 
arc interested in what is essentially the simplest model 
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of prcdator-prcy interactions, wc include only birth pro- 
cesses BE BB, death processes A ^ E, B ^ E, 
and predator-prey interactions AB ^ AA, AB ^ AE. 
Here {b,di,d2,Pi,P2) are rate constants. The symbol E 
corresponds to what would be available sites in a spatial 
model. In this non-spatial model, the i?'s are {N — n—m) 
passive constituents of the system, which are required for 
prey reproduction, and which result in intra-specific prey 
competition. Note, the overall number of A, B and E 
constituents is fixed to be N. The dynamics of the model 
consists of choosing constituents at random and imple- 
menting the rules given above. The time dynamics of the 
model can either be numerically simulated or studied an- 
alytically using the formalism of master equations Q ■ 
In the latter case the transition rates T(n', m'ln, m) from 
the state (n, m) to the state (n', m!) are given by 



PLM, takes the form 
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where the b and di have been scaled by a factor of {N— 1) 
and the di by a factor of N. 

We stress that the individuals of a given species in our 
model are identical, and thus the term ILM should not 
be confused with "agent based models" which are often 
designed to study the ecological effects of behavioral and 
physiological variation among individuals. We have al- 
ready given an extensive discussion of this approach else- 
where in the context of competition models 1^, and we 
refer the reader to this paper for a fuller discussion of the 
formalism. For predator-prey models however, including 
stochastic dynamics leads to more marked effects than in 
competition models, as we now discuss. 

The master equation for the probability that the sys- 
tem consists of n predators and to prey at time t, 
P(n, TO, t), is 



dP{n, TO, t) 
dt 
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where the step operators £ are defined by their actions 
on functions of n and m by f{n, m, t) = /(nil, to, t) 
and £y^f{n, m, t) ^ f{n, to ± 1, t). 

The mean field limit of this ILM may be obtained by 
multiplying |(2J) by n and to in turn, and subsequently 
summing over all allowed values of m and n. This gives 
equations for the mean values /i = {n)/N and /2 = 
(to) /N in the limit — > cx) if we ignore terms which are 
down on others and make the replacements (to^) — > 
(to) 2 and (mn) — > {m){n). This mean field theory, or 
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The eqs. © are frequently referred to as the Volterra 
equations, to distinguish them from the Lotka- Volterra 
equations which have no term in /2//!r Q . The constants 
/i, r, and K are simply functions of the rate constants: 



fi = di , r = 2b — d2 
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and the linear numerical and functional responses are 
given by n(/2) = 2pi/2 and g(/2) = 2{'pi + P2 + 6)/2 
respectively. 

As is well known , the analysis of this model shows a 
complete absence of cycles. There is a single fixed point 
for which the predators and prey have non-zero popula- 
tion sizes. Denoting these stationary values by ^\ and 

is] 

/2 , then in terms of the original rate constants they are 
given by: 
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The stability of this fixed point may be studied may per- 
forming linear stability analysis. This results in a stabil- 
ity matrix which is given by 
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FIG. 1: Predator and prey densities as a function of time. 
The upper panel shows the predator density /i for A'^ = 3200. 
The blue line is calculated from numerical integration of the 
mean field Volterra Eqs. Q. The purple line is the average of 
the predator density time series from 500 replicates generated 
from the ILM, and is almost indistinguishable from the mean 
field solution. The red line is the predator density time series 
for a single typical replicate. The lower panel is the equivalent 
plot for the prey density /2. Parameter values are b = 0.5, 
di = 0.1, d2 = 0.0, pi = 0.5, and p2 = 0.1. 
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We have expressed the entries in terms of the fixed point 
values, since these are manifestly positive, and it is easy 
to see that the eigenvalues of A both have a negative real 
part, implying that the fixed point is stable. The entries 
of A will appear again below in the analysis of the cy- 
cling behavior for finite N. For now let us simply remark 
that, while there is no limit cycle in the Volterra system 
(O, a limit cycle does exist in the Lotka- Volterra equa- 
tions (obtained by taking K — > oo), but it is neutrally 
stable due to a conserved quantity in the model. This 
unrealistic behavior disappears with the introduction of 
a finite carrying capacity, K, in but, as mentioned 
above, leads to a complete absence of cycling behavior. 
We will show below that cycles can be found in the ILM, 
but only when N is finite; the TV — > oo limit which was 
taken in order to derive the PLM, eliminates the cycles 
present in the original ILM. 

To see this, let us first note that the ensemble aver- 
aged population density of the ILM, determined from 
numerical simulations, agrees beautifully with the solu- 
tion of this deterministic model (Fig. 1, purple and blue 
lines respectively) showing a decaying oscillatory tran- 
sient followed by a constant steady-state density, typical 
of a Volterra system. In marked contrast, individual real- 
izations of the ILM show large persistent cycles (Fig. 1, 
red line). The amplitude of these oscillations is much 
larger than the naive estimate based on the law of large 
numbers. In fact, the oscillations are of order (I/a/ZV) 
as would be expected, but amplified by a very large fac- 
tor due to a noise-induced resonance effect, as explained 
below. 

This cycling behavior can be investigated analytically 
by extracting an "effective theory" valid for large N, 
through applying a standard method to the master equa- 
tions, due to van Kampen 01 . Essentially the method in- 
volves the replacements n/N = fi + x/\/N and m/N = 
f'l + y/V^ in the transition probabilities that appear 
in the master equation. By changing from a description 
based on the (discrete) variables n and m to one based 
on the (continuous) variables x and y, terms of different 
orders in can be identified in the master equation: 
the leading order terms gives rise to a deterministic set of 
equations and the next-to-leading order terms give rise to 
a linear Fokker-Planck equation. The leading order set of 
equations (mean field theory) are the PLM and are the 
Volterra equations J^l which we have already obtained 
by a more direct method. 

At next-to-leading order, rather than write down the 
Fokker-Planck equation, it is simpler to give the set of 
Langevin equations to which it is equivalent 0- They 
take the form 
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These are a pair of differential equations which describe 
the stochastic behavior of the ILM at large N: x{t) and 
y{t) are stochastic corrections to the deterministic be- 
havior of the predator and prey densities respectively. 



at large but finite N. The constants, aij, appearing in 
Eq.© are exactly the entries of the matrix A, Eq. © 
found from a linear stability analysis about the non- 
trivial fixed point of Eq. Q . The noise covariance matrix 
hij, which is responsible for generating the large-scale 
oscillations, cannot be determined from Eq.® and is de- 
rived from the master equation using the van Kampen 
expansion. Since the noise is white, bij = {f}i[ijj)f)j{-'Uj)) 
is independent of the frequency lo. The explicit expres- 
sions for these constants are 
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As discussed earlier, it is not the average behavior of 
replicates that interests us, but rather measures which 
characterize the oscillations. Examining x and y as func- 
tions of frequency allow us to determine the nature of the 
oscillations. 

To search for oscillations in noisy data, one of the most 
useful diagnostic tools is the power spectrum P{lu) = 
(|i(w)p), where x{uj) is the Fourier transform of x{t). 
Taking the Fourier transform of Eq. lO, solving for x{lu), 
and averaging its squared modulus, we find 
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where a and /3 are functions of the ILM rates: a = 
foiia|2 + 2fei2ai2|a22| +&22I12 ^nd (3 = bn. The constants 
in the denominator have the especially simple forms: 
^0 ~ 1^12 10-21 1 a-nd r = 1 022 1- The spectrum predicted 
by Eq. © gives the blue line shown in Fig. 2. The agree- 
ment with the spectrum obtained from simulation of the 
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FIG. 2: A plot of the power spectrum P{uj) for the preda- 
tor time series, as a function of frequency uj. The red line 
corresponds to P calculated from 500 replicate runs of the 
ILM. The blue line is the prediction from our theory, namely 
Eq. Q . The parameter values are the same as those described 
in the caption to Fig. 1. The inset shows the analogous power 
spectra (data and theory) for the prey time series. 
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ILM (red line) is excellent. Note, the naive 0{l/y/N) es- 
timate of the size of stochastic fluctuations corresponds 
to the zero frequency value of P{l<j). Fig. 2 clearly illus- 
trates the very large amplification of these fluctuations 
due to the resonance effect. 

The spectrum given above is reminiscent of that for a 
simple mechanical system — namely a linear damped 
harmonic oscillator, with natural frequency f2o a-nd 
driven at frequency oj. In a mechanical oscillator the 
driving frequency must be tuned to achieve resonance. 
In the stochastic predator-prey model described here no 
tuning is necessary. The system is driven by white noise, 
as shown in Eq. ((TJ, which covers all frequencies — thus 
the resonant frequency of the system is excited without 
tuning. We stress, the noise which drives the system is 
not external, but arises from the demographic stochas- 
ticity contained in the individual processes which define 
the model. We also stress that there is no external or 
environmental stochasticity in our model, and that the 
resonance phenomenon we report here is not related to 
"stochastic resonance" . 

The damping term, represented by the constant F, lim- 
its the amplitude of the oscillations. Predator-prey sys- 
tems for which F happens to be small will be at risk of 
extinction through resonant oscillations, despite having 
large population sizes. The resonant oscillation occurs in 
the regime 2ai2|a2i| > 032, where the resonant frequency 
tJo = -^/Og — F2/2 is real. A similar analysis can be car- 
ried out to obtain the spectrum for the prey time series. 
It again has the form Eq. 0, but now with a = buo^i 
and (3 = 622. The positions of the peaks for these two 
power spectra are only weakly dependent on the a's and 
/?'s, and so they are almost coincident. 

Predator-prey systems (and related host-pathogen sys- 
tems) have been studied theoretically for decades. Most 
of the previous studies have focused on the role of envi- 
ronmental stochasticity, the relevance of non-linear inter- 
actions or of spatial effects, to explain the mechanism of 
cycling d i. Ill IH IH El lli . Some authors have dis- 
cussed the role that demographic stochasticity may have 
on cycles 0, . Most of this discussion has been quali- 



tative; the nearest to our own discussion was a prescient 
analysis by Bartlett nearly fifty years ago in which 
he postulated equations similar to 0. However, he did 
not note the existence of a resonance, and so proposed 
cycles with an amplitude which were not enhanced by 
this effect, and which were therefore of limited biological 
interest. 

The idea that external perturbations with a dominant 
frequency can entrain the predator-prey dynamics in a 
cyclic nature is fairly intuitive; the phenomenon we dis- 
cuss here is more fundamental and less intuitive. The 
noise in our system is internal — purely a result of the de- 
mographic stochasticity inherent in discrete birth, death, 
and prcdation events. The key point is that this internal 
noise has a flat power spectrum in the frequency domain; 
in other words, it excites all frequencies of the system si- 
multaneously. These excitations are typically of limited 
interest in a large population of N individuals, since they 
give rise to small 0(1 /^/N) fluctuations about the mean 
population densities. Such is the case, for example, in 
competition models @|. The predator- prey system, and 
related ones such as epidemic models, are exceptional, in 
that the equations describing linear fluctuations about 
the steady-state are susceptible to resonant amplifica- 
tion in the vicinity of an internal frequency fl, which 
is a property of the population itself. The internal noise, 
in exciting all frequencies, automatically resonates the 
system giving rise to large oscillations in the population 
densities. This phenomenon is "emergent" in the truest 
sense of the word. We expect that this resonance mecha- 
nism will occur in other stochastic systems in which the 
mean field theory shows damped oscillations. 
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